HPAVC/HPAVC CD-ROM.iso/pc/EARTHQUA.ZIP/EQSIM/EQSIM.CPP

C/C++ Source or Header | 1995-02-04 | 19.9 KB | 901 lines

// EQSIM.CPP - The simulator for the earthquake damage prevention program.
// Generates the building movement and displays it.
// Created by Misha Koshelev.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#include <dos.h>
#include <string.h>
#include <alloc.h>
#include <graphics.h>
#include "svga256.h"
#include "gendisp.h"

// Misc. definitions for init_globals
#define WAVELET_NUM 15
#define DELTA_OMEGA 0.2
#define FRACTAL_DIM 0.35

// Definitions for active and hybrid control
#define DELAY 3
#define CONTROL_NORMAL 2.0
#define HYBRID_SPOINT 1.0

// Program defaults
#define FELASTIC 1
#define FFRICTION 1
#define TSTEP 0.1
#define SNUM 300
#define BMASS 10.00
#define MAGNIFICATION 1.0

// Some defines
#define EXIT 255.255

typedef double (* cntrlfunc)(double);            // The definition for the
                                               // control function.
typedef struct
{
    double f, a, v, x, e;                    // The force, accelaration,
                        // veloctity, energy used,
                        // and x movement of a
                        // building
                        // where the original position
                        // is 0.0
} bms;

bms far *blmov = NULL;                          // The building movement
                                               // array.
double felastic_cos, ffriction_cos;             // The coefficients for
                                               // elastic and friction
                                               // forces.
double tstep;                                                       // The time step used for
                                               // the building movement.
int   snum;                                                                    // The number of steps.
double bmass;                                                               // The mass of the building.

// Predefined headers
int huge Svga256Detect(void);
double no_cntrl(double);
double active_cntrl(double);
double hybrid_cntrl(double);
double execute_cntrl(cntrlfunc, double);
void init_qglobals(double, double, double, double, int, double);
void init_gen_qglobals(double, double, double, int, double);
void done_qglobals(void);
void done_gen_qglobals(void);
void gen_blmov(cntrlfunc);
void plot_blmov(int, int, int, int, int, int, int);
void plot_disp(int, int, int, int, int, int, int);
void loadpal(char *);
void hrainbow(int, int, int, int, int, int, int, int);
void vrainbow(int, int, int, int, int, int, int, int);
double total_blmov(void);
double max_blmov(void);
double total_eused(void);
double show_control(cntrlfunc, char *, double);
void show_control_txt(cntrlfunc, char *, double);
double select_intensity(void);

// Global variables
int bpoly[10];

// Fill pattern, to fill the building with
char bpattern[8] = {0xFF, 0xFF, 195, 195, 195, 195, 0xFF, 0xFF};

// Set if program should write the result to a file (default)
int wtof = TRUE;

FILE far *tmf, far *mmf, far *ef;

// Set if demo mode (default off)
int demomode = FALSE;

// Set if no graphics (default off)
int graphics = TRUE;

// The main program.
//
#ifdef STANDALONE
void main(int argc, char *argv[])
#else
void eqsimmain(void)
#endif
{
   double ti, i, intensity;
   int ai;
   int argprocessed = FALSE;
   int gd = DETECT, gm, errorcode;

#ifdef STANDALONE
   // Check for parameters
   if (argc > 0)
   {
      for (ai = 1; ai < argc; ai++)
      {
     // If parameter is disable write to file, disable it
     if (!strcmp(argv[ai], "/NOOUTPUTFILES"))
     {
        argprocessed = TRUE;
        wtof = FALSE;
     }
     // If parameter is demo mode, enable it
     if (!strcmp(argv[ai], "/DEMO"))
     {
        argprocessed = TRUE;
        demomode = TRUE;
     }
     // If parameter is no graphics, disable it
     if (!strcmp(argv[ai], "/NOGRAPHICS"))
     {
        argprocessed = TRUE;
        wtof = TRUE;
        demomode = FALSE;
        graphics = FALSE;
     }
     if (argprocessed == FALSE)
     {
        printf("Error: Invalid parameter %s", argv[ai]);
        exit(1);
     }
     argprocessed = FALSE;
      }
   }
#endif

#ifdef STANDALONE
   if (graphics)
   {
      // Install the Super VGA 256 color driver
      installuserdriver("Svga256", Svga256Detect);

      // read the result of installation
      errorcode = graphresult();

      if (errorcode != grOk)  // an error occurred
      {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     printf("Press any key to quit...");
     getch();
     exit(1);
      }

      // Initalize graphics mode
      initgraph(&gd, &gm, "");

      // read result of initialization
      errorcode = graphresult();

      if (errorcode != grOk)  // an error occurred
      {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     printf("Press any key to quit...");
     getch();
     exit(1);
      }

      // Load the pallette
      loadpal("std.pal");
   }
#else
   if (graphics)
      cleardevice();
#endif

   // If we should write to files, open the files
   if (wtof)
   {
      tmf = fopen("totalmov.eqs", "wt");
      mmf = fopen("maxmov.eqs", "wt");
      ef = fopen("energy.eqs", "wt");
   }

   if (wtof && (!tmf || !mmf || !ef))
   {
      printf("Error: Could not create data files");
      exit(1);
   }

   // Set up the TEX header
   if (wtof)
   {
     fprintf(tmf,"\\settabs 4 \\columns\n\n");
     fprintf(tmf,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
     fprintf(mmf,"\\settabs 4 \\columns\n\n");
     fprintf(mmf,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
     fprintf(ef, "\\settabs 4 \\columns\n\n");
      fprintf(ef,"\\+Intensity&No Control&Active Control&Hybrid Control\\cr\n");
   }

   if (graphics)
   {
      // Initialize the polygon's coordinates
      bpoly[0] = 300;
      bpoly[1] = 249;

      bpoly[2] = 340;
      bpoly[3] = 249;

      bpoly[5] = 200;
      bpoly[7] = 200;

      bpoly[8] = bpoly[0];
      bpoly[9] = bpoly[1];
   }

   for (i=TSTEP; i<9.5; i+=TSTEP)
   {
NewIntens:
      // Write the intensity to the data files:

     if (wtof)
     {
        fprintf(tmf, "\\+%.1f&", i);
        fprintf(mmf, "\\+%.1f&", i);
        fprintf(ef, "\\+%.1f&", i);
     }

      // No Control
      intensity = pow(10, i/3);
      init_qglobals(intensity, FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
      if (!graphics)
      {
         show_control_txt(no_cntrl, "No Control:",i);
      }
      if (graphics)
      {
        ti = show_control(no_cntrl, "No Control:",i);
        if (ti == EXIT)
           goto Done;
        if (ti != 0.0)
        {
           done_qglobals();
           i = ti;
           loadpal("std.pal");
           goto NewIntens;
        }
      }
      done_gen_qglobals();

      // Show Active Control
      init_gen_qglobals(FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
      if (!graphics)
      {
         show_control_txt(active_cntrl, "Active Control:",i);
      }
      if (graphics)
      {
        ti = show_control(active_cntrl, "Active Control:",i);
        if (ti == EXIT)
           goto Done;
        if (ti != 0.0)
        {
           done_qglobals();
           i = ti;
           loadpal("std.pal");
           goto NewIntens;
        }
      }
      done_gen_qglobals();

      // Show Hybrid Control
      init_gen_qglobals(FELASTIC, FFRICTION, TSTEP, SNUM, BMASS);
      if (!graphics)
      {
         show_control_txt(hybrid_cntrl, "Hybrid Control:",i);
      }
      if (graphics)
      {
        ti = show_control(hybrid_cntrl, "Hybrid Control:",i);
        if (ti == EXIT)
           goto Done;
        if (ti != 0.0)
        {
           done_qglobals();
           i = ti;
           loadpal("std.pal");
           goto NewIntens;
        }
      }
      if (wtof)
      {
        fprintf(tmf, "\\cr\n");
        fprintf(mmf, "\\cr\n");
        fprintf(ef, "\\cr\n");
      }
      done_qglobals();
   }

   if (wtof)
   {
      fclose(tmf);
      fclose(mmf);
      fclose(ef);
   }
Done:
#ifdef STANDALONE
   closegraph();
#endif
}

#ifdef STANDALONE
// Returns what video mode to use.
//
int huge Svga256Detect(void)
{
   return SVGA640x480x256;
}
#endif

// The function for no control
//
// Input: Time
//
// Output: 0
//
double no_cntrl(double t)
{
   return 0.0;
}

// The active control function.
//
// Input: Time
//
// Output: The number to add to the force
//
// The algorithm used is:
//
// f        = -hx(t - T)
//  control
//
// where h is a coefficient, t is the current time, and T depends on the
// delay of the control system.
//
double active_cntrl(double t)
{
   double r = (-CONTROL_NORMAL * blmov[(int)(t / tstep) - DELAY].x);

   blmov[(int)(t / tstep)].e = r * (blmov[(int)(t / tstep)].x - blmov[(int)
                 (t / tstep) - 1].x);
   return r;
}

// The hybrid control function.
//
// Input: Time
//
// Output: The number to add to the force
//
// The algorithm used is:
//
// if x(t - T) <= p then:
//
// f        = 0.0
//  control
//
// else
//
// f        = -hx(t)
//  control
//
// where h is a coefficient, t is the current time,
// and T depends on the delay of the control system.
//
double hybrid_cntrl(double t)
{
   double r = 0.0;

   blmov[(int)(t/tstep)].e = 0.0;
   if (blmov[(int)(t / tstep) - DELAY].x > HYBRID_SPOINT)
      r = (-CONTROL_NORMAL * blmov[(int)(t / tstep)].x);
   return r;
}

// Executes the control functions if the array element is greater than
// DELAY-1.
//
double execute_cntrl(cntrlfunc cf, double param)
{
   if ((int)(param / tstep) - DELAY < 0)
   {
      blmov[(int)(param/tstep)].e = 0.0;
      return 0.0;
   }
   return cf(param);
}

// Initializes the earthquake parameters.
//
// Input: Intensity, Elastic force coefficient, Friction force coefficient,
// the time step, the number of steps, and the building's mass.
//
// Output: None
//
void init_qglobals(double intens, double fe_cos, double ff_cos, double t_step,
     int s_num, double b_mass)
{
   snum = s_num;
   if (blmov)
      done_qglobals();
   blmov = (bms far *)farcalloc(snum, sizeof(bms));
   if (!blmov)
   {
      closegraph();
      printf("Not enough memory\n");
      exit(1);
   }
   for (int j=0; j<snum; j++)
   {
      blmov[j].x = 0.0; blmov[j].v = 0.0; blmov[j].f = 0.0; blmov[j].a = 0.0;
      blmov[j].e = 0.0;
   }
   tstep = t_step;
   felastic_cos = fe_cos;
   ffriction_cos = ff_cos;
   bmass = b_mass;
   init_globals(intens, WAVELET_NUM, tstep * snum, DELTA_OMEGA, FRACTAL_DIM);
}

// Initializes the building movement parameters without initializing the
// earthquake parameters.
//
// Input: Elastic force coefficient, Friction force coefficient,
// the time step, the number of steps, and the building's mass.
//
// Output: None
//
void init_gen_qglobals(double fe_cos, double ff_cos, double t_step,
     int s_num, double b_mass)
{
   snum = s_num;
   if (blmov)
      done_gen_qglobals();
   blmov = (bms far *)farcalloc(snum, sizeof(bms));
   if (!blmov)
   {
      closegraph();
      printf("Not enough memory\n");
      exit(1);
   }
   for (int j=0; j<snum; j++)
   {
      blmov[j].x = 0.0; blmov[j].v = 0.0; blmov[j].f = 0.0; blmov[j].a = 0.0;
      blmov[j].e = 0.0;
   }
   tstep = t_step;
   felastic_cos = fe_cos;
   ffriction_cos = ff_cos;
   bmass = b_mass;
}

// Destroys the global variables
//
void done_qglobals(void)
{
   farfree(blmov);
   blmov = NULL;
   done_globals();
}

// Destroys the global variables without destroying the earthquake parameters.
//
void done_gen_qglobals(void)
{
   farfree(blmov);
   blmov = NULL;
}

// Generates the building movement
//
// Input: A control function
//
// Output: None
//
// The following algorithm is used for every step where t is the current time.
//
// x(t) = v(t - s) * s + x(t - s)
//
// v(t) = a(t - s) * s + v(t - s)
//
// f(t) = f         + f        + f         + f
//         external    elastic    friction    control
//
// f        = -kx(t)
//  elastic
//
// f         = -nv(t)
//  friction
//
// a(t) = f(t)
//        ----
//         m
//
// e(t) = f        * x(t) - x(t - s)
//         control
//  /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\
// for active control and 0 for all others.
//
// where s is the time step, m is the mass of the building, k and n are
// coefficients, f(t) is the force in a current time, x(t) is the position
// related to the original position at the current time, v(t) is the velocity
// at the current time, and a(t) is the accelaration at the current time.
//
void gen_blmov(cntrlfunc cf)
{
   // Generate values for time 0
   blmov[0].x = 0.0;
   blmov[0].v = 0.0;
   blmov[0].f = disp(0.0) + execute_cntrl(cf, 0.0);
   blmov[0].a = blmov[0].f / bmass;
   // Generate all of the other values
   for (int j=1; j<snum; j++)
   {
      blmov[j].x = blmov[j-1].v * tstep + blmov[j-1].x;
      blmov[j].v = blmov[j-1].a * tstep + blmov[j-1].v;
      blmov[j].f = disp(j * tstep) + (-felastic_cos * blmov[j-1].x) +
                      (-ffriction_cos * blmov[j-1].v) +
                      execute_cntrl(cf, j * tstep);
      blmov[j].a = blmov[j].f / bmass;
    }
}

// Plots the building movement.
//
void plot_blmov(int x, int y, int x1, int y1, int borcolor, int normcolor,
           int limitcolor)
{
   int cx, cy, middle, color;
   int limit;

   setcolor(borcolor);
   rectangle(x, y, x1, y1);       // Draw border rectangle
   middle = ((y1-y) / 2) + y;
   limit = x1+1;
   if ((x1+1)-x > snum)
      limit = x+snum;
   for (cx=x;cx<limit;cx++)
   {
      cy = (int)blmov[cx-x].x * MAGNIFICATION + middle;
      color = normcolor;
      if (cy < y)
      {
     cy = y;
     color = limitcolor;
      }
      if (cy > y1)
      {
     cy = y1;
     color = limitcolor;
      }
      putpixel(cx, cy, color);
   }
}

// Plots the earthquake displacements
//
void plot_disp(int x, int y, int x1, int y1, int borcolor, int normcolor,
           int limitcolor)
{
   int cx, cy, middle, color;
   int limit;

   setcolor(borcolor);
   rectangle(x, y, x1, y1);       // Draw border rectangle
   middle = ((y1-y) / 2) + y;
   limit = x1+1;
   if (x1-x > snum)
      limit = x+snum;
   for (cx=x;cx<limit;cx++)
   {
      cy = (int)disp((cx - x) * tstep) + middle;
      color = normcolor;
      if (cy < y)
      {
     cy = y;
     color = limitcolor;
      }
      if (cy > y1)
      {
     cy = y1;
     color = limitcolor;
      }
      putpixel(cx, cy, color);
   }
}

#ifdef STANDALONE
// Loads the palette from a raw palette file.
//
void loadpal(char *fn)
{
   int color;
   int r, g, b;
   FILE *pal;

   pal = fopen(fn, "rb");

   if (pal == NULL)
   {
      closegraph();
      printf("Error: Could not open palette file %s", fn);
      exit(1);
   }
   for (color = 0; color < 256; color++)
   {
      fread(&r, 1, 1, pal);
      fread(&g, 1, 1, pal);
      fread(&b, 1, 1, pal);
      setrgbpalette(color, r, g, b);
   }
   fclose(pal);
};
#endif

// Calculates the total building movements
//
double total_blmov(void)
{
   double sum = 0.0;

   for (int j=0; j<snum; j++)
      sum += fabs(blmov[j].x);

   return sum;
}

// Returns the maximum building movement.
//
double max_blmov(void)
{
   double max;

   if (snum == 1)
     return fabs((double)blmov[0].x);
   max = (fabs(blmov[0].x) > fabs(blmov[1].x)) ?
                         fabs(blmov[0].x)
                          : fabs(blmov[1].x);
   if (snum == 2)
      return max;
   for (int j=2; j<snum; j++)
     max = (max > fabs(blmov[j].x)) ? max : fabs(blmov[j].x);
   return max;
}

// Returns the total energy used.
//
double total_eused(void)
{
   double sum = 0.0;

   for (int j=0; j<snum; j++)
      sum += fabs(blmov[j].e);

   return sum;
}

// Shows the control wave and the earthquake wave.
// Returns 0.0 or to what number to change the current intensity.
//
double show_control(cntrlfunc cf, char *cfname, double intens)
{
   char tmp[20];
   char ch;
   int y;

   // Set the colors
   setfillstyle(SOLID_FILL, BLUE);
   bar(1,1,640,300);

   setfillstyle(SOLID_FILL, 79);
   bar(1,250,640,300);

   // If in demo mode, write the word DEMO
   if (demomode)
   {
      sprintf(tmp, "DEMO");
      settextstyle(TRIPLEX_FONT, HORIZ_DIR, 8);
      setcolor(32);
      outtextxy(300, 10, tmp);
      setcolor(WHITE);
      settextstyle(DEFAULT_FONT, HORIZ_DIR, 1);
   }

   // Show the intensity
   sprintf(tmp, "Intensity: %.1lf", intens);
     outtextxy(10, 10, tmp);
   y = 10 + textheight(tmp) + 2;

   // Show the earthquake displacement
   sprintf(tmp, "Earthquake Disp:");
     outtextxy(130, 300, tmp);
   plot_disp(10,310,310,479,WHITE,WHITE,RED);

   // Show the building movement
   gen_blmov(cf);
   outtextxy(440, 300, cfname);

   plot_blmov(320,310,620,479,WHITE,WHITE,RED);

   // Display the data
   sprintf(tmp, "Total: %.3lf", total_blmov());
     outtextxy(10,y,tmp);

   if (wtof)
      fprintf(tmf, "%.3lf&", total_blmov());

   y += textheight(tmp);
   sprintf(tmp, "Max. Move.: %.3lf", max_blmov());
     outtextxy(10,y,tmp);

   if (wtof)
      fprintf(mmf, "%.3lf&", max_blmov());

   y += textheight(tmp);
   sprintf(tmp, "E. Used: %.3lf", total_eused());
     outtextxy(10,y,tmp);

   if (wtof)
      fprintf(ef, "%.3lf&", total_eused());

   y += textheight(tmp);

   bpoly[4] = 340;
   bpoly[6] = 300;

   setfillstyle(SOLID_FILL, 128);
   fillpoly(5, bpoly);
   setfillpattern(bpattern, WHITE);
   fillpoly(5, bpoly);

   setfillstyle(SOLID_FILL, BLUE);

   for (int k=0; k<SNUM; k++)
   {
      if (bpoly[4] != (int)blmov[k].x + 340 && bpoly[6] !=
       (int)blmov[k].x + 300)
      {
     bar(1,y+1,640, 249);
     bpoly[4] = (int)blmov[k].x + 340;
     bpoly[6] = (int)blmov[k].x + 300;

     setfillstyle(SOLID_FILL, 128);
     fillpoly(5, bpoly);
     setfillpattern(bpattern, WHITE);
     fillpoly(5, bpoly);
     setfillstyle(SOLID_FILL, 124);

     setfillstyle(SOLID_FILL, BLUE);
      }
      else
       delay(17);      // If no time taken to draw delay for 17 ms
   }

   if (!demomode)
   {
#ifdef STANDALONE
      outtextxy(200, 265, "Press any key to continue...");
      outtextxy(200, 275, "To change intensities press <I>");
      outtextxy(200, 285, "To exit press <ESC>");
#else
      outtextxy(200, 265, "Press any key to continue...");
      outtextxy(200, 275, "To exit press <ESC>");
#endif
      while (!kbhit());
      ch = getch();
      setfillstyle(SOLID_FILL, BLACK);
      bar(1,1,640,480);
#ifdef STANDALONE
      if (ch == 'i' || ch == 'I')
     return select_intensity();
#endif
      if (ch == 27)
      {
     done_qglobals();
     if (wtof)
     {
        fclose(tmf);
        fclose(mmf);
        fclose(ef);
     }
#ifdef STANDALONE
     closegraph();
#endif
     return EXIT;
       }
   }

   if (demomode)
   {
      if (kbhit())
      {
     done_qglobals();
     if (wtof)
     {
        fclose(tmf);
        fclose(mmf);
        fclose(ef);
     }
#ifdef STANDALONE
     closegraph();
#endif
     return EXIT;
      }
      setfillstyle(SOLID_FILL, BLACK);
      bar(1,1,640,480);
   }
   return 0.0;
}

// Shows the current intensity, etc. and writes the data to a file.
//
void show_control_txt(cntrlfunc cf, char *cfname, double intens)
{
   char tmp[20];
   char ch;

   // Show the intensity
   sprintf(tmp, "Intensity: %.1lf", intens);
   printf("%s", tmp);
   printf("\n");

   printf("%s\n", cfname);

   cprintf("Please wait...");

   // Generate the building movement
   gen_blmov(cf);

   // Display the data
   sprintf(tmp, "Total: %.3lf", total_blmov());
   printf("\n%s\n", tmp);

   if (wtof)
      fprintf(tmf, "%.3lf&", total_blmov());

   sprintf(tmp, "Max. Move.: %.3lf", max_blmov());
   printf("%s\n", tmp);

   if (wtof)
      fprintf(mmf, "%.3lf&", max_blmov());

   sprintf(tmp, "E. Used: %.3lf", total_eused());
   printf("%s\n\n", tmp);

   if (wtof)
      fprintf(ef, "%.3lf&", total_eused());

   if (kbhit())
   {
      ch = getch();
      done_qglobals();
      if (wtof)
      {
     fclose(tmf);
     fclose(mmf);
     fclose(ef);
      }
      exit(0);
   }
}

#ifdef STANDALONE
// Asks the user for a new intensity.
//
double select_intensity(void)
{
   double r;
   int gd = DETECT, gm;

   restorecrtmode();
   printf("What intensity do you want to switch to (X.X)? ");
   scanf("%3lf", &r);
   if (r > 9.5)
     r = 9.5;
   setgraphmode(getgraphmode());
   return (r);
}
#endif